function A = jacobian_f_robot(x0, u, dt)

eps = 1e-3;

for j=1:numel(x0)
   x = x0;
   x(j)=x0(j)+eps;
   y2=f_robot(x, u, dt);
   x(j)=x0(j)-eps;
   y1=f_robot(x, u, dt);
   
   A(:,j)=(y2-y1)/(2*eps);
end
